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Abstract 

We investigate the thermodynamics and dynamics of the electroweak phase tran- 
sition by modelhng the infrared physics with classical Yang-Mills Higgs theory. We 
discuss the accuracy of this approach and conclude that, for quantities whose deter- 
mination is dominated by the infrared, the classical method should be correct up to 
parametrically suppressed (ie 0{a)) corrections. For a Higgs self-coupling which at 
tree level corresponds to mu — 50GeV, we determine the jump in the order parameter 
to be = 1.5gT, the surface tension to be a = 0.07g'^T^, and the friction coefficient 
on the moving bubble wall due to infrared bosons to he rj = P/v^ = 0.03 ± .OOAg^T'^. 
We also investigate the response of Chern-Simons number to a spatially uniform 
chemical potential and find that it falls off a short distance inside the bubble wall, 
both in equilibrium and below the equilibrium temperature. 
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1 Introduction 



Since Sakharov proposed that the universe's baryon number could have been gener- 
ated dynamically early in the Big Bang many different particle physics scenar- 
ios have been proposed to realize this mechanism. One of the most interesting is 
based on the observation made 20 years ago by t'Hooft that baryon number is vi- 
olated already in the minimal standard model [0; and by the later realization that 
the violation is very efficient at temperatures above the electroweak phase transition 
temperature, where the Higgs condensate is absent @, §]. It was subsequently shown 
that electroweak baryogenesis could actually occur in minimal extensions of the stan- 
dard model involving one or more additional Higgs doublets, and that the resulting 
asymmetry was naturally of the order of magnitude required in cosmology P|. 
Other, potentially more efficient baryogenesis mechanisms were proposed 0, |^ , both 
in these theories and in other extensions of the standard model. What is most ex- 
citing about these scenarios is that they link the origin of matter in the universe to 
physics beyond the standard model that could be probed by the next generation of 
Higgs-seeking accelerators, and other experiments looking for novel sources of C and 
CP violation. 

Because of these developments there has been quite intense interest in the elec- 
troweak phase transition. But the quest to understand baryogenesis at the electroweak 
scale has been bedeviled by difficulties. It is now believed that we know what ba- 
sic issues are involved. As temperature drops, the Higgs field develops a condensate 
abruptly in a first order phase transition. This proceeds by the nucleation and sub- 
sequent growth of bubbles of the low temperature, broken symmetry phase. Outside 
the bubbles, baryon number is violated, but inside it is conserved. The plasma de- 
parts from equilibrium on and very near the boundaries of these bubbles, and if there 
is some CP violation, for instance in the Higgs sector, then all necessary conditions 
for the creation of a baryon number excess exist; and it would be preserved to this 
day because the rate of baryon number violation is very low inside the bubbles, if the 
phase transition is suitably strong. It is also strongly believed that the communica- 
tion of the C and CP violating physics to the infrared gauge-Higgs fields, which are 
responsible for the baryon number violation, should be conducted by the fermions, 
and that the transport of the fermions in the presence of the wall may be relevant. 
But turning this understanding into predictions requires an accurate description of 
the dynamics of the phase transition, which has so far been poorly developed. The 
dynamics, in turn, cannot be well understood until the thermodynamics are under 
control; and all aspects require a complete knowledge of the underlying electroweak 
physics, which of course is still lacking. 

In the coming years we will hopefully learn the required information about elec- 
troweak physics, for instance the exact nature of the Higgs mechanism, the existence 
or absence of low energy supersymmetry, and details of any non-CKM CP violation. 
But our ignorance of these details does not prevent work on the thermodynamics and 
dynamics of the electroweak phase transition from proceeding. The tools developed to 
investigate the simplest case, the minimal standard model, should be straightforward 
to extend to more complicated models; in fact the extension may often be analytically 
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tractable, as seems to be the case generically for the thermodynamics 0, p!0[| . 

One tool which has proven less useful than hoped is perturbation theory. While 
one loop perturbation theory can determine the phase transition temperature with 
0{a) accuracy |]1T|, it is less reliable for the details of the phase transition, such as the 
jump in the order parameter, the latent heat, the surface tension, and so forth. This 
became clear when Arnold and Espinosa computed the effective potential to two loops 
||T^ , 13|. At one loop a term of form —g^cfy^T/n arises which generates a Higgs vev 
squared of 0^ ~ gf^T^/A^Tr^; and at two loops a new term of form —g'^cfP'T'^ ln(0/T)/7r^ 
arises, which alone would generate 0^ ~ (yf^T^/Avr^. As far as the strength of the 
phase transition is concerned, perturbation theory is at best an expansion in A/^f^, or 
{mH/mwY. 

Because A receives large radiative corrections, it is difficult to make the ratio \/ 
very small; we also know that, if the minimal standard model is correct, the ratio 
cannot be small, because direct searches for the Higgs boson have ruled it out below 
a mass of 60GeV. It therefore seems reasonable to take X ^ g^ parametrically. This is 
also natural from the point of view of unit analysis, if we allow h to be dimensionful. 
To do otherwise would require a rearrangement of perturbation theory away from 
the usual loopwise expansion, which can be understood in the vacuum theory as an 
expansion in powers of h. We see from the above discussion, though, that at the 
phase transition temperature, for quantities such as the jump in the order parameter, 
the perturbation series is not an expansion in h, which helps to explain its poor 
convergence. 

Farakos et. al. have elucidated the reason for this, by showing how the infrared 
thermodynamics of the minimal standard model at finite temperature is very well 
described by a "dimensionally reduced" 3 dimensional field theory [Q. this analysis 
has allowed an efficient nonperturbative investigation of the phase transition which 
is now bearing fruit [^, 0. 

An at first sight unrelated development is the idea that the rate of anomalous 
baryon number violation in the high temperature phase of the classical theory is 
computable numerically, and that the rate in the classical theory should be the same, 
up to 0{a) corrections, as the rate in the quantum theory [|r^, |r^, |TU|. This has 



allowed the first quantitative measurement of a dynamical, infrared property of finite 
temperature Yang-Mills Higgs theory above the phase transition temperature. The 
validity of this idea is fortified by the observation that the thermodynamics of the 
classical theory coincides with the thermodynamics of the quantum theory in the 
approximation of dimensional reduction ^ . 

In fact we believe that the success of the dimensional reduction technique, and 
the nonperturbative nature of the infrared quantum theory, arise precisely because 
the infrared bosonic modes attain large occupation numbers which cause the theory 
to mimic a classical theory. We will argue in Section |^ that this mimicry extends 
to the dynamics, so that all infrared dominated dynamical properties of the theory 
which possess cutoff independent limits in the classical theory will take the same 
value, up to 0{a) corrections, in the quantum theory. We then review an evolution 
algorithm for the lattice cut off, classical field theory in Section |^. Since the classical 
theory has the same thermodynamics as the quantum theory in the dimensional 
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reduction approximation, this evolution algorithm is also a microcanonical Monte- 
Carlo algorithm for the dimensionally reduced theory; so we develop and apply the 
tools for using it to investigate the thermodynamics of the phase transition in Section 
^. Then we turn to dynamics; in Section ^ we use the classical technique to compute 
the contribution of infrared bosons to the friction felt by a moving bubble wall as it 
sweeps through the plasma, in a near equilibrium approximation. In Section ^ we 
investigate baryon number violation in the presence of an out of equilibrium bubble 
wall. Section |^ concludes. There are also two appendicies. Appendix A presents the 
tools for the numerical study of Chern-Simons number (Ncs) motion, and investigates 
baryon number violation in each phase; in particular it presents evidence that the 
observed rate of baryon number violation in the broken phase arises from ultraviolet 
lattice artefacts. Appendix B presents the details of how to extract the surface tension 
of the bubble wall. 

We make no serious attempt to extrapolate to the continuum limit either for the 
thermodynamical or dynamical properties, so the results presented here should be 
considered rough and preliminary; the emphasis is on development of techniques. 



2 Discussion of the classical approximation 

As discussed in the introduction, we will investigate the use of classical Yang-Mills 
Higgs theory as a surrogate for the Standard Model. In important respects the clas- 
sical theory does not resemble the quantum theory at all (for instance, it has an 
infinite, or cutoff dependent, heat capacity); it should only be used for those dy- 
namical and thermodynamical properties which are dominated by infrared physics, 
where it should give results which reproduce the parametrically leading term in the 
full quantum theory. Hence, for instance, the classical theory has been used to inves- 
tigate baryon number violation in the symmetric electroweak phase, a phenomenon 
thought to be infrared dominated |TP|, |2T|, It has also been used, 

in a slight disguise, to investigate the thermodynamics of the phase transition. 

To see the latter point, consider the dimensional reduction program of Farakos 



et al |2J, |T5|, ^ 0, a systematic, semiperturbative approach to determining the 
strength and other thermodynamic properties of the phase transition. Their idea 
is the following. All thermodynamic properties can be derived in the Matusbara 
formulism, ie by considering a Euclidian path integral in which time runs from to 
j3 = 1/T with periodic boundary conditions for bosons and antiperiodic boundary 
conditions for fermions. To compute the thermal expectation value of an operator (9, 
we find 

^ ' /P$I?A/^exp(-5) ' ^ ' 

S = j dxo J (fx{C + Cct) , 
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The subscript 4 means that these values depend on renormahzation point exactly as 
they do in the 4 dimensional vacuum theory, as do the wave functions; the theory 
also requires counterterms represented by £ct- (For simplicity we have not written 
terms involving fermions, but they should be present.) By Fourier transforming the 
time direction, one obtains a 3 dimensional theory in which the Fourier components 
with nonzero Matsubara frequency appear as a Kaluza-Klein tower of massive modes. 



Following |25|, they construct an infrared effective theory of the zero frequency mode, 
integrating out the massive modes by computing a set of correlators in each theory 
and matching them in the infrared. Provided that the operator O consists of spatially 
extended, equal time combinations of bosonic operators, its expectation value is then 
well approximated by 

^ JV^VA^Oexpi-(3H) 
^ ' JV^VA>^exp{-(3H) ' ^ ' 

1 1 2 

^A2 + m^o$t$ + A($t$)2. (4) 

Here the wave functions and couplings do not renormalize but have their values fixed 
(at a given temperature) by the matching process. (We have dropped an extremely 
small Aq term and a slight correction to the coefficient in the y4Q$^$ term, as in ||15|| .) 



The mass terms, on the other hand, do depend on the renormahzation procedure, 
as indicated by the subscript 0, meaning the bare values. They are related to the 
renormalized values, determined in the perturbative matching procedure, by 

mlifi) = m^Q + 6ml{fi) , (5) 

and similarly for mjj. In regulations where linear divergences do not vanish, such 
as lattice regulations, the counterterm is substantial and positive, so the bare mass 
squared may need to be small or negative. 

We have deliberately used different notation than Farakos et. al. to emphasize 
that the path integral is over an action which looks like H/T, with H "almost" the 
Hamiltonian of the classical theory. 

To investigate the relationship between the dimensionally reduced theory above 
and the thermodynamics of the classical bosonic theory, we follow a line of reason- 
ing developed in ||20[. In the classical theory, thermodynamics are described by the 



partition function 

Z = j VA1VE^VYW^5 [{DiEif + g'^Reli^^^^ exp{-f3H) (6) 

H = lQi^^.i^^. + iEfEf)+ntn+(A$)^A$+ (7) 



^2 3$"^$ + A(<I)"^$)^ 

where the delta function enforces Gauss' law. It can be written by introducing an 
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integration over a Lagrange multiplier Aq and adding to the Hamiltonian 

The measure for and 11 is now trivial and the integrals are Gaussian. Performing 
them, the partition function becomes 

Z = J V^VA^exp{-(3H) (9) 

+ m^o^t* + A($t$)2 , (10) 

identical to Eq. except that m^yQ is forced to be zero. The actual Debye mass 
squared then equals the counterterm. In lattice regulation, the counterterm turns out 
to be [13 

S = 3.17591, (11) 

A-Ka 

which grows linearly with 1/a. One choice is to set the lattice spacing a so that 
this is actually the correct Debye mass, but this is not essential. All that is really 
needed is that the screening be efficient enough that the influence of the Aq field on 
the infrared physics should be perturbatively computable (which is the case for the 
physical value of m^), and is true for the value in Eq. ( pA]) for reasonable values of a). 
In this case we can use the results of to relate the thermodynamics of the lattice 
system to the thermodynamics without the Aq field, which in turn can be related to 
the thermodynamics with the Aq field and the appropriate Debye mass. Except for 
the need for this correction, an investigation of the thermodjTiamics of the classical 
theory is equivalent to an investigation of the dimensionally reduced theory. 

The preceding discussion suggests a profound connection between the thermody- 
namics of the classical theory and the nonperturbative infrared physics of the quantum 
theory. To explore whether the connection extends to real time properties, we will 
briefly investigate the real time perturbative expansion of each theory. 

Let us start with classical field theory. For simplicity we shall take as usual the 
example of scalar field theory, and we assume some regularization, like lattice 
regularization, is present. The classical thermal ensemble is defined as: 

where O is some quantity of interest. The classical field 0j(x) and its momentum 
7ri(x) are those at some particular initial time tf. the measure is to be thought of 
as a measure on the space of initial conditions for the field. The Hamiltonian H 
is Hq + Hint, with Hq = J cPx^{7rf + (pf) and Hint = ^ / d^xXcjyf. Since both the 
measure and the Hamiltonian are time independent, the expectation value (0)/3 is 
independent of ti. As far as O is concerned, it may be any function of the classical 
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field and momentum, evaluated at any time. For example O = 0(x, t)0(O,O) gives 
the classical unequal time two point correlator, where 0(x, t) is the classical solution 
defined by the initial conditions 0j(x) and 7ri(x) at t = tj. 

Solving the classical theory perturbatively is straightforward. If one is interested 
in the correlator (0(x, t)0(O, 0)), for example, one can choose ti = and solve the 
classical field equation 

(d^ - V' + m')^ = -^A0' (13) 

with an integral equation: 

Cj){x, t) = ^freeix, ^) " ^ dt' J £x'Gr{x - x\t - t')(j)%x' , t') . (14) 

where 0/ree is a solution of the free theory, 

G„(. 1') = I ,15) 

J (27r)'^ ujk 



is the retarded Greens function, and uJk = \^W+m?. The iteration of this equation 
produces the solution for (p{x,t) to all orders in A. Of course (j)freeix,t) is easily 
expressed in terms of 0(x, 0) and 7r(x, 0): 

j^e^^'"" (0(k, O)cos(t<;fct) + 7r(k, 0)sm{ujkt)/ujk) . (16) 

To evaluate (0) one expands e~^^'"* in powers of A, and performs the Gaussian 
integrations over 0(x, 0) and 7r(x, 0) occurring at each order in A. These integrals are 
summarized by the generating functional: 

^Jd''x^f,.e.{x)J{x)s^^^^^^ = ^\Sd^xjd^x'.]{x)G(x,x')J(x) ^-^^^ 
G{x,x')^<<Pfree{x)<Pfree{x')> = T / -^^^''^^"'^'^ • i^^) 

To summarize, the interactions occur in two places: first, in e~^^ , which defines 
the thermal state and determines the equal time correlators, and second, in the clas- 
sical evolution of the fields between the times of interest in unequal time correlators. 
The same is true in the quantum theory. 

In the quantum theory, one wants to compute the quantum expectation values of 
operators: 

iOh - (19) 

where the trace is over any complete set of states, and O is the product of field 
operators, expressed in terms of Heisenberg fields. As in the classical theory, one can 
evaluate this expression in an interaction picture defined at any time (e.g. t = 0). 
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Rewriting (^) in terms of the interaction picture operators (f)\'x,t) and U{t) 
Te~ ^ ^0 '^^'"^«t (*') we have 



Tr ( e-'3^-(0^(x, t)0^(O, 0) + 0^(0, 0)0^ (x, t)) 



[i\U{-ihf3)U\t)(f)\x,t)U{t)(p\0,0)\t) + c. c. (20) 



where the sum is over a complete set of interaction picture occupation number states. 
This expression is straightforward to evaluate, by expanding all the U operators in 
powers of A and then evaluating the ensuing free field correlators. 

To see the connection with the classical theory, note that W{t)(f)\:K,t)U{t) is just 
the original Heisenberg field, which obeys the classical field equation (|13|), and as in 
the classical theory, this may be solved by iterating the integral equation (0), with 
the identical retarded Greens function. The perturbative expansion of the Heisenberg 
operator yields this Greens function via [0^(x, t), 0^(x, t')] = ihG r^x — x' , t — t') , t > t'; 
the h cancels the in the time evolution operator, to all orders in A. So the 
'dynamical' part of the calculation of unequal time correlators in the classical and 
quantum theories are actually identical^. It is also remarkable that when organized 
this way, both theories have exactly the same set of Feynman diagrams, it is just the 
Feynman rules that are different. 

Differences arise in two places. First, in the expectation values of the resulting 
series of free fields. In the quantum theory the generating function is given as in (16) 
but with Gciass{x,x') replaced by 



G. 



quantum\-^ 1 



(2^ 



3ik.(x— x') 



COs(co'fct) 



+ 



(21) 



(This 'thermal Wicks theorem' has been discussed recently by |^ - the simplest way 
to derive it is to note that the relevant path integral is Gaussian, from which (16) 
follows. The 'thermal Wicks theorem' is true for any initial density matrix which 

is Gaussian.). The second difference arises because 

not exactly equal to e~/?ifj„t(0). This difference can be attributed to the difference 
between the quantum and classical (or 'dimensionally reduced') thermal states. 

We now see when the classical and quantum results agree. For low frequency 
bosonic modes at high temperatures huo/T -C 1, the occupation number is high and 
the bracketed expression in Eq. ( PT| ) can be expanded in an asymptotic series. 



1 

T+2 



T 



+ 



1 huok 



(22) 



•^except that one must be careful about operator ordering in using the Heisenberg operator equa- 
tion of motion. The products of operators one gets wiU not be completely ordering averaged. 
However, re-arranging their order only introduces commutators which are 0(h), which by dimen- 
sion counting means 0{huj/T). It turns out that when one asks questions about ordering averaged 
products of operators, as we do here, these re-ordering corrections actually first appear at second 
order in h, see Kql. 
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and the leading term precisely reproduces the classical Greens function. Note that 
the usual h in the propagator is cancelled by a in the occupation number. Second, 
as noted above, H-^^{t) is not equal to H-^^{0) - the imaginary time dependence of 
the free field operators causes a difference in the thermal state corrections. However 
H-^^{t) can be expanded as a Taylor series in r, and the r dependent corrections 
come in the form hukT < huj^/T . Thus again for huj/T <^ 1, the corrections are small 

B 

To summarize, as long as one is interested in phenomena which are dominated by 
low frequency bosonic modes, for which huj/T <^ 1, the the classical field theory pro- 
vides a good description. Note that the thermal classical theory contains much more 
than just tree diagrams - it sums up all loop diagrams as well, in the approximations 
first that the quantum propagator is replaced by the classical {% independent) piece 
and second that the thermal state is taken to be that of the 'dimensionally reduced' 
three dimensional theory. 

The above argument works perfectly in a regularized (e.g. lattice) field theory. The 
only catch in applying it to a continuum field theory is that the condition huj/T 
cannot be true for the very high k modes. For large k the quantum and classical 
propagators are very different, and the ultraviolet behavior of the theory approaches 
that of the vacuum quantum theory. However, if we are interested in long wavelength 
correlation functions, these ultraviolet modes will only emerge in diagrams in which 
there are a few closely spaced vertex insertions, which we should be able to replace 
with an expansion in local operators. The classical theory will only make sense as a 
regulated, infrared effective theory with a Lagrangian which takes into account these 
local, quantum effects. For instance, the couplings will be renormalized, to a scale 
given roughly by the temperature. 

In addition to short distance vacuum corrections, which simply lead to the usual 
renormalization of the couplings to a scale given roughly by T, there are thermal (non- 
vacuum) short range corrections, whose form is not restricted by Lorentz invariance 
because the thermal bath chooses a preferred time direction. There has been intense 
research into how these effects influence the behavior of the infrared (soft) modes, 
and it has been shown that the parametrically leading effects can be summarized as 
a set of "hard thermal loop" effects which can be incorporated into the Lagrangian 
of the infrared theory . 



For this reason, Bodeker et. al. proposed to include the hard thermal loop effects 
in numerical investigations of the classical field theory |2^. In fact, their work shows 



that they are already included, because the most ultraviolet classical excitations in 
a (lattice cut off) classical simulation perform the same role as the high frequency 
thermal excitations in the quantum theory; all that is different from the quantum 
theory is the shape of the cutoff and the total strength of the hard thermal loop effects. 
The total strength of the hard thermal loop effects generally cancels in calculations 
of dynamical properties of the plasma |^ , because the hard thermal loops represent 



"^In fact for equal time correlators of only, these thermal state corrections occur only at order 
(huJk/T)^ , as is seen by writing U[—ihfi) ~ g+PHo^-PH -y^ji^j-^ (quantum) corrections arising from 
commutators of Hq with Hint. At first order in h there is only the single commutator, which is odd 
in TT and therefore does not contribute to the correlator. 
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both a bath of particles against which to scatter, and a bath of particles which screen 
interactions so that only shorter range scatterings can take place. It is not generally 
believed that the total number of particles available to contribute to the hard thermal 
loops affects the infrared dynamics, and if it does, then dynamical properties would 
show a power law in a (the lattice spacing) dependence. Testing for a small a limit 
to dynamical properties therefore constitutes a check that the magnitude of hard 
thermal loops does not matter, or at least that the dynamical property in question 
converges to a limit in the (parametrically justified) limit of large hard thermal loop 
contributions. The only remaining concern is that the functional form of the hard 
thermal loops is incorrect because of cutoff artefacts; but by varying the form of the 
cutoff one can test for this as well, and at least in the case of the motion of Chern- 
Simons number it appears that lattice artefacts are small and consistent with zero 

Hence, we conclude that the classical theory can be used to examine the thermo- 
dynamics of Yang-Mills Higgs theory and allows the study of the infrared dynamics 
as well, though here there are some legitimate concerns involving hard thermal loops. 



3 Numerical implementation of the theory 

The numerical implementation of classical, real time Yang-Mills Higgs theory is well 
developed in the literature |^T], |2D|, The implementation we use here is identical 
to that of we consider the theory at zero Weinberg angle (no f/(l)y) and take 
as degrees of freedom an SU(2) matrix on every link of a 3+1 dimensional lattice 
which is toroidal in the 3 space directions but infinite in the time direction, and a 
fundamental complex scalar at each vertex. The lattice action is 



-yi--Trf/n + — ^— y 1 - -Trf/n, 

tr 2 ^ (At)2tr 2 



-EE ^(*(^' t) - U^{x)^x + I, t)y{<l>{x, t) - U,{x)<l>{x + z, t)) 

x,t i 

+ ^ E ^(*(^' t) - t^o(x)$(x, t + At))t($(x, t) - Uo{x)^x, t + At)) 



2 \ 



x,t 



(23) 



We absorb the gauge coupling into the lattice temperature Pl, which is related to the 
lattice spacing and the continuum temperature through 

Pl = ^. (24) 

We also give the Higgs fields the same wave function normalization as the gauge 
fields, which is natural and computationally convenient. The Higgs field is treated 
as four independent real entries, and the relation between the lattice value and the 
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continuum one is 

/?l$t<|. = Accent = ^-^^^cont (25) 

9 9 

for the normal continuum definitions of (f) and $. Witli tliis wave function normal- 
ization, the scalar self-coupling is related to the usual continuum one by 

Xl = ^^/9\ (26) 

which is parametrically order 1 if A ~ (7^, as is natural from the renormalization 
structure of the theory. The only small numbers in the classical theory are oc a 
the lattice spacing, and the renormalized Higgs mass squared. 

Varying the action with respect to the link in the time direction generates a con- 
straint, analogous to the continuum condition DiEI/g"^ = Re(n^ir"$) = p, which is 
Gauss's law; variation with respect to the space links and the Higgs field components 
give equations of motion which allow all future times to be determined from two 
neighboring time slices (provided that these initial conditions satisfy the constraints), 
once an ambiguity in the time evolution, due to the freedom to change gauge inde- 
pendently at each spacetime point, has been used up by choosing the gauge which 
always makes the time links Uq = I the identity. 

In practice we keep track of the values of fields on one time slice, and "momenta" 
U{t-At/2) = $(t)-$(t-At) and Ef{t-At/2) = -l/2Trir'^f/i(t)f//(t - At) (after 
setting the Uq to J); as long as E is small, the latter relation can easily be inverted 
to update U. 



We thermalize the system with the algorithm developed in [^. That is, beginning 
from an arbitrary initial condition, we repeatedly draw the momenta from the correct 
thermal distribution (achieved by choosing them as Gaussian random variables, and 
then orthogonally projecting to the constraint surface) and then evolve the system 
under the equations of motion for some time, allowing the thermalization to mix with 
the coordinates U and $. The momenta are then discarded and the procedure is 
repeated, as many times as desired. The algorithm has a time stepsize ambiguity. 



which we handle as in |2^; the thermalization is only accurate to 0((At)^), a level 
which is sufficient because the evolution algorithm is also inaccurate at this level. In 
this work we always use At = 0.05 in lattice units, which is sufficient to hold stepsize 
systematics below statistical errors. 

As discussed above, a thermalization algorithm for classical Yang-Mills Higgs the- 
ory can be considered a canonical ensemble (fixed temperature) Monte-Carlo algo- 
rithm for the dimensionally reduced theory. (In fact, except for the Gauss constraint, 
it exactly resembles the hybrid algorithm of Euclidean lattice gauge theory.) The al- 
gorithm is very efficient at exploring the thermodynamics of one phase; but it is very 
bad at going between the two phases. This is because it must move smoothly from 
one phase to the other, and mixed phase configurations contain phase boundaries 
which have positive surface tensions. In the limit of large box size, the suppression 
of such configurations grows roughly as exp(— 2L1L2/0"), with Li and L2 the two 
shortest lengths of the box and a the surface tension. It is necessary to change phase 
several times to get good statistics on the free energy difference of the two phases. 



11 



and hence to determine the critical temperature, but as the box size grows, it will 
become essentially impossible for the canonical evolution algorithm to do so; it will 
only be capable of thoroughly exploring one minimum, but not of comparing the two. 
For this reason, the literature generally considers a multicanonical ensemble (in which 
a global reweighting term is added to make mixed phase configurations more favor- 
able but is then accounted for in computing thermal averages of operators) a more 
powerful technique for exploring the thermodynamics of first order phase transitions. 

In fact, properties of the metastable phases, the equilibrium temperature, and 
virtually all other themodynamic properties can be extracted by using a microcanon- 
ical (fixed energy) ensemble. A microcanonical ensemble is achieved by thermalizing 
the system to some temperature, but then allowing it to evolve under the equations 
of motion indefinitely, without ever re-randomizing the momenta. The system is 
strongly ergodic, so for general initial conditions it will thoroughly explore the fixed 
energy subsurface of phase space. In the large volume limit, this would become ap- 
proximately equivalent to a canonical ensemble, except that there is a first order 
phase transition. There is a finite range of energies where the fixed energy equilib- 
rium configuration is mixed phase, and if one can find one mixed phase configuration 
in this range, then the microcanonical algorithm can use it to thoroughly explore 
mixed phase configurations, and in particular to extract the phase transition temper- 
ature and information about the phase interfaces. Since the whole range of mixed 
phase configurations occur at a single temperature, the canonical ensemble is not well 
suited to exploring phase coexistence. This difference between the two ensembles is 
illustrated in Figure ||. 

To exploit this property of the evolution algorithm, it is necessary to find a way to 
measure temperature during a fixed energy evolution, and to very gradually increase 
or decrease the energy so that a range of configurations can be explored. If there were 
no Gauss constraint, it would be easy to measure the temperature. The momenta 
would be true Gaussian degrees of freedom, and their average energy (averaged over 
the set of momenta and over time to remove fluctuations) should obey equipartition. 
The Gauss constraint only mildly complicates this picture. The constraints are linear 
in the momenta, and remove one Gaussian degree of freedom each. Although they are 
not local, so we do not know an orthogonal basis for the remaining set of independent 
Gaussian degrees of freedom, we still know their total number, so the total kinetic 
energy should equal (number of degrees of freedom=10)T/2. Averaging over time 
removes statistical fluctuations and gives a clean value for the system temperature. 

It is also quite easy to gradually change the amount of energy in the system. Our 
algorithm to do this is to simulate a gradual, adiabatic expansion or contraction of 
the lattice. All coordinates {U and $) are updated as usual, but all momenta E, U 
are multiplied by the same factor (1 + eAt) each time step. The system heats with 
time constant 1/e. If e is negative, the system will cool. Since the Gauss constraint 
is linear in E and 11, this heating algorithm identically preserves Gauss' Law. 

It is also sometimes necessary to make the heating or cooling more local, so that 
the temperature can be kept uniform even if some phenomenon liberates heat locally. 
To do this we bin the lattice into boxes or slabs, and measure the temperature in 
each. Then each box has its momenta multiplied using a different e, chosen to drive 
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each box independently towards some desired temperature. The algorithm generates 
small violations of Gauss' Law on the interfaces between boxes, which we remove with 
the orthogonal projection algorithm presented in |^ . 

Note that neither of these heating algorithms will keep the system in equilibrium; 
of course no heating algorithm will. But if they are applied gradually, the system 
should remain very close to equilibrium (though it will tend to superheat or supercool 
into metastable phases, as noted earlier). 



4 Thermodynamic properties 

We are now ready to explore the equilibrium properties of Yang- Mills Higgs theory. 
We will not attempt a thorough investigation of the strength of the phase transition 
as a function of A^, as very accurate calculations already exist ||16|, |3^; rather we will 
investigate how well the microcanonical technique can be applied to determining the 
various properties of the phase transition. Similarly, we have not attempted to make 
a small lattice spacing extrapolation; all the data presented below are for /S^, — 8 or 
6 and Xl = 0.20, which, in the notation of 0, is x = 0.05. 



4.1 Metastability and hysteresis 

The first thing we can investigate is the temperature dependence of order parameters, 
including metastable branches. We do this by thermalizing a 30'^ lattice with bare 
Higgs mass squared mjjQ = —0.3223 in lattice units and Xl = 0.20 at a temperature 
above the phase transition temperature, estimated to occur where m'jjf) equals the 
one loop counterterm, 

_ (9 + 6A.)S 

which happens when (3l ~ 8.0. The system is gradually cooled, and the lattice 
temperature and order parameters such as are averaged in time bins longer 
than the lattice length but much shorter than the full length of the cooling. After 
the order parameter jumps to the broken phase value, the system is heated back to 
the original temperature. As a function of temperature, order parameters exhibit 
hysteresis, returning to the symmetric phase at a higher temperature than they left 
it. This is the signature of a first order phase transition. The hysteresis curve for 
<l'^$, a once smoothed defined in Appendix B, and the traces of a few sizes of 
Wilson loops are presented in Figures ^ and ^. Smoothing $ greatly reduces the 
contributions from the most ultraviolet fluctuations, but barely touches the infrared 
fluctuations. We expect that almost all of the value of in the symmetric phase 
should arise from ultraviolet fluctuations, which should contribute Pl^/tt ~ 8.1 to 
fl^, so the value of should be much lower in the symmetric phase after 
smoothing; but almost all of the difference between phases should be infrared and 
unaffected by smoothing. Figure |^ verifies this, and also shows that almost all of the 
random thermal fluctuations in the order parameter are infrared, because the details 
of the fluctuations in are almost unchanged by the smoothing. Note also that. 
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except near the spinodal point, the fluctuations in the symmetric phase value of 
are much smaller than in the broken phase; this is because of interference between 
fluctuations and the condensate in the broken phase. Also note that there is a large 
correlation between the fluctuations of different sized Wilson loops, and between the 
fluctuations in Wilson loops and the fluctuations in $^<l>. 

Examining the Wilson loop plots, we see that, in the symmetric phase, traversing 
a 9 X 9 loop yields an almost completely random SU(2) phase, so the symmetric 
phase is disordered on the scale of Pl. In the broken phase traces of Wilson loops 
fall off more slowly, but 13 x 13 Wilson loops show almost no order and it cannot 
be meaningful to speak of anything in a nontrivial representation of SU(2) as having 
any correlations beyond this scale. Hence the box size was abundantly larger than 
the longest possible correlation length and the results should represent the inflnite 
volume limit. We also verifled this by repeating the run on a 20^ lattice; the results 
were the same within error, but the fluctuations were larger because they were not 
averaged over as much four- volume, so we will not present the results here. 

4.2 Equilibrium temperature 

The results presented in Figure ^ can be used to determine the jump in $t$ at any 
temperature for which both phases are reasonably metastable. However, it cannot be 
used to determine the nucleation rate, because in the early universe the supercool- 
ing occurred much more gradually than in any concievable simulation and the true 
nucleation point occurs when the tunneling probability is too small to observe in a 
simulation. It also cannot give us the equilibrium temperature, although if we knew 
the equilibrium temperature we could use the hysteresis plot to flnd the jump in the 
order parameter at equilibrium. To get the equilibrium temperature, we must estab- 
lish phase coexistence in a microcanonical evolution and measure the temperature. 

We do this as follows. We thermalize an x x 192 rectangular box with 
symmetric boundary conditions at the same value of and rrij^Q used above and 
(3l = 8.0. Then we apply a small perturbation to wi^q, with an amplitude which varies 
sinusoidally along the long direction of the box. The symmetric phase is favored in one 
region and the broken phase is favored in another, so a mixed phase conflguration 
is established. The perturbation is then slowly removed, and as it is removed the 
system is heated or cooled by a thermostat which tries to balance / ^^^d^x at a value 
intermediate between the two phases. Once the perturbation is completely lifted, the 
system is allowed to evolve for a long time (at least 500 lattice lengths) without any 
heating or cooling to equilibrate fully in the mixed phase conflguration and erase all 
record of the process by which a mixed conflguration was generated. It is then evolved 
for a long period of time, again without any heating or cooling, during which the 
energy in kinetic degrees of freedom is averaged and used to establish the equilibrium 
temperature. We have checked that the box we used was abundantly longer than 
that required to contain two domain walls and a region of each phase, and that the 
system remained in the two phase conflguration during the whole run. The two phase 
nature can be seen clearly by averaging optionally applying several iterations 
of smoothing, over the two short directions of the lattice and plotting against the 
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long direction, as shown in Figure |. The reason that fluctuations in phase boundary 
positions do not drive the system to one or the other phase is that a fluctuation 
which expands the broken phase hberates latent heat, which raises the temperature 
and makes the symmetric phase more favorable, and similarly a fluctuation which 
expands the symmetric phase absorbs latent heat and makes the broken phase more 
favorable. 

We measured the equilibrium temperature in a 16 x 16 x 192 box and in a 32 x 
32 X 192 box; the answers are within errorbars and average to I3l = 8.059 ± .0020. 
The jump in the order parameter at this temperature, read off the hysteresis 

plot, is 9.5 ± .2. This corresponds to a jump in the continuum order parameter 
of 1.54gT {g the weak coupling). Two loop perturbation theory predicts a jump 
of l-24g, which is smaller^; this might at flrst seem surprising, considering that the 
values found by Kajantie et. al. [|T^] for A^, on either side of = 0.2 are closer to 
the two loop perturbative value. The reason is that they perform an extrapolation to 
zero a (inflnite Pl) based on data at several values of (3l, whereas the result we quote 
above is for one value of (3^ and contains flnite lattice spacing artefacts. Because 
Ai: is small and the jump in the order parameter is quite sensitive to its value, the 
most important of these effects may be those which shift the effective value of A^,. 
One arises because we have Aq flelds with flnite Debye mass squared; if we integrated 
them out, then they would shift A^ by 

-Sg^T 



SX = 5Xl = — , „ ~ -0.0188J8//3l . (28) 

There is also a linear in a correction to Xl arising from one loop diagrams. The 
one loop contribution to the effective potential is 

Vi{(f)o)=T. mlim)dm, (29) 

^ Jm{<f>=0) 

where the sum is over massive degrees of freedom and lim), the one loop 3d lattice 
tadpole graph, is computed in [ll5| . 
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+ 0(a^), ^ = 0.1529. (30) 



The leading, 1/a term generates the linearly divergent mass squared correction, the 
next term gives the negative cubic term which determines the order of the phase 
transition, and the arn^ term produces an 0(a) correction to A which, summed over 
degrees of freedom, is 



^In the thermalization algorithm, E(t) is drawn from the Gaussian distribution and £'(t + At/2), 
the value required for the leapfrog algorithm, is determined by a half leapfrog step; but the tem- 
perature quoted is the sum of -^/l — (At)^i?^(t + At/2)/ At, so the value we quote will have a weak 
stepsize dependence. 

^Here and throughout we use Eq. (34) of |]l5| , with all terms involving dropped, for the two 
loop effective potential 
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which equals —0.014 for (3l = 8. Hence, the simulations described correspond to of 
roughly 0.167 rather than 0.2, and the two loop perturbative estimate for is 1.40^', 
quite close to the actual value. 

Note that accounting for finite a shifts in as above does not completely remove 
finite a or even all linear in a errors, because we have not corrected wave functions 
or removed high dimension operators. It is probably not profitable to pursue high 
statistics calculations until the required corrections have been computed. 

Also note that the jump in together with the value for m^Q, determines 

the latent heat, as shown in [0. A more direct measure of the latent heat is the 
temperature change during the spinodal jump from one phase to the other, which can 
be looked up on the hysteresis plot (though it must be remembered that the heating 
or cooling of the system continued during the jump). The value obtained from the 
jump in temperature is consistent with the latent heat obtained from the jump in 

4.3 Surface tension 

Another interesting property of the equilibrium system is the surface tension of the 
phase interface. The surface tension can be determined from the power spectrum of 
fluctuations of the interfaces. To understand this, consider a surface with a very large 
surface tension; since it is taut, it should be flat. But for flnite surface tension, the 
entropy associated with having nonzero fluctuations on the surface prevents it from 
being perfectly flat. For long wavelength fluctuations, the fluctuation amplitude is 
small compared to the wavelength, and different fluctuations approximately decouple; 
on average they are populated according to equipartition. Thus the infrared limit of 
the power spectrum (square of the Fourier coefficients) of the bubble surface deter- 
mines the surface tension. The details are given in Appendix B, where we also discuss 
how we define the bubble wall surface. 

To apply the technique discussed in the appendix, we evolved an equilibrated 
mixed phase configuration for on order 1000 lattice lengths, recording both bubble 
wall surfaces every 2 lattice lengths of time. Each surface is Fourier transformed, and 
the square of each Fourier coefficient is averaged over the run. For each value of 
there are several independent coefficients (from the real and imaginary parts of one 
or more Fourier coefficients from two walls), and we average these and take the error 
bars to be the standard deviation of the mean. 

We plot the resulting power spectrum, multiplied by in? the square of the Fourier 
mode number, for a 32 x 32 x 192 lattice at /5l ~ 8 and for a 36 x 36 x 144 lattice at 
/3l ~ 6 in Figure ^ This is not a log-log plot; the departure from behavior in 
going from = 1 to = 25 is only on order a factor of two. Part of this departure 
from strict power law behavior may arise from the smoothing involved in defining 
the bubble wall surface, and some of it may represent interactions between the high 
frequency modes on the wall. To make the infinite wavelength extrapolation we fit 
the data with an exponential; the fit has x^/v of 0.67 and 1.4 for the 32^ and 36^ 
cross section cases and yields surface tensions in physical units of .0681 ± .QQQQg^T^ 
and .0739 ± .00085''^T^ respectively. Note that for /^^ = 6 the finite a systematics are 
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different; the earlier estimate for the correction to Xl gives = 
one should have expected a larger surface tension. 

The one loop analytic estimate of the surface tension is 



.160 in this case, so 



2^1 (0), (32) 

where Vi is the one loop effective potential. For Xl = 0.167 the result is .0254(yf^T^ 
and for = .160 it is .0285'^T^. One may also apply the same formula using the 
two loop effective potential, though this is slightly inconsistent since one is still using 
the tree level kinetic term without one loop wave function corrections. The results 
are .OSOg'^T^ and .0865'^T^ respectively. The large difference between the one and 
two loop perturbative estimates reflects the 0^ dependence of the surface tension. 
Since the two loop perturbative value for (p and the lattice value are quite close, it 
is not surprising that the two loop surface tension is quite close to the lattice value. 
The lattice is larger than the two loop value, but the lattice surface tension is 



smaller; this may be the beginning of a trend of low surface tensions found in p6 
though it is difficult to say until the remaining 0{a) effects are accounted for. A lower 
surface tension is also the direction one would expect from including wave function 
corrections, as discussed in [p^ . 



We conclude that the microcanonical technique is an efficient and promising way 
of extracting all interesting thermodynamical properties of the phase transition; it 
should be further pursued after a more careful accounting of 0{a) corrections has 
been made. 



5 Friction on the bubble wall 

The previous section merely uses the Hamiltonian evolution of classical Yang-Mills 
Higgs theory as a microcanonical Monte-Carlo algorithm for thermodynamic investi- 
gation, but as stressed in Section |^, the classical Hamiltonian evolution should also 
give information about dynamics, provided that the physics involved is infrared dom- 
inated. This includes two previously elusive phenomena, friction on the bubble wall 
from infrared bosons and the motion of Chern-Simons number near a moving bubble 
wall. 



5.1 General discussion 

The velocity attained by a moving bubble wall during the cosmological electroweak 
phase transition is one of the key ingredients for models and calculations of baryon 
number production. Several authors have considered the problem [R3|, R3, 153, 



39| , ^ p] , 1^, and quite a bit is known. Two effects prevent the bubble wall 
from "running away" and establish its terminal velocity; frictive effects arising from 
the departure from equilibrium of massive species due to the motion of the wall, and 
hydrodynamic effects arising from the liberation of latent heat. The hydrodynamic 
properties of the plasma are dominated by thermal energy particles, which hold almost 
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all of the energy and momentum of the plasma, so these effects cannot be studied by 
classical techniques. Fortunately, except perhaps for energy transfer across the bubble 
wall by ballistic leptons, the hydrodynamics are well under control |3^, ^ . 

The friction from massive particles depends on a high power of the particle mass, 
and so the only important particles for consideration are the top quarks and the bosons 
of the Yang-Mills Higgs system. Friction from top quarks cannot be calculated in the 
classical theory, but it arises at a higher parametric order than the friction from 
bosons, where the leading parametric contribution is infrared dominated and should 
be reproduced in the classical simulations. 

To see this, we begin by stating what we mean by the friction. We will only be 
interested in this paper with the friction in the case that the departure from equilib- 
rium is small, so that a fluctuation dissipation formula can relate it to measurable 
equilibrium correlators In this case we define the friction coefficient as 



ri = lim — , 



(33) 



where v^, is the average velocity of a planar bubble wall when the pressure difference 
between phases is P. As we will discuss below, this friction coefficient is related to the 
diffusion constant for the random motion of the equilibrium interface, and converting 
from lattice units to thermal units establishes the parametric behavior of the classical 
1] as r]ci oc g^T"^. 

The friction coefficient for the quantum theory can be computed at lowest order 
in a loopwise expansion from a fluctuation-dissipation theorem and the answer 
is equivalent to the friction from free scattering, thermal particles on the bubble wall 
PP| , which has been calculated in 
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where the first integral arises from particles which scatter from the wall and the second 
is from particles which fly over the wall from each side. If we make the approximation 
l/(exp(/3i?) — 1) ~ {1/ [3E), which is precisely the classical approximation to the Bose- 
Einstein population factor, then the integrals give precisely the flrst 0{g^(jr') term. 
Hence we can understand the result as a classical part, which (recalling that ~ gT) 
is order g^T'^, plus a correction which is 0{g^\n.l/ g). The most important point is 
that the result is totally dominated by the infrared contribution, and extending the 
classical approximation for the Bose distribution out to arbitrarily large momenta 
only makes a parametrically suppressed error; hence the friction from bosons is (at 
leading parametric order) an infrared, classical effect. This means, of course, that 
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the loopwise expansion will be unreliable, because the infrared is strongly coupled; 
instead we should use classical real time simulations to determine the friction. 

For the fermions, Eq. ( ^If ) applies but with the Bose-Einstein statistics replaced 
by Fermi-Dirac statistics l/(exp(/3i?) + 1) (and of course the 6 replaced by a 12). The 
resulting friction is 



r/ ~ 12 



(36) 



1287r2 \4 C/y0 

which begins at 0{g^lnl/g). Friction from top quarks is parametrically suppressed 
and will not appear in the classical theory, which is obvious because the classical 
theory is purely bosonic. The top quark contribution to friction may still be signif- 
icant, simply because it depends on gy and gy is numerically much larger than g. 
The absence of a large infrared contribution should also make the top quark friction 
computable, although the naive perturbative series must be resummed into Boltz- 



mann equations to account for on-shell near singularities and may then require 
simplifying analytical approximations to make the calculation feasible |^T], 

Fortunately, most of the top quark friction arises from the small departure from 
equilibrium of thermal energy particles; because Fermi statistics are well behaved in 
the infrared, the departure from equilibrium of top quarks will have a small (paramet- 
rically suppressed) influence on the infrared bosonic sector. The two frictions should 
be additive with 0{a) error, and an investigation of the friction in the bosonic theory 
is well motivated. 



5.2 Fluctuation dissipation relation 

As mentioned above, the 0{g^) contribution to the friction coefficient can be com- 
puted in the limit of a small departure from equilibrium by a fluctuation dissipation 
argument, as follows. Consider an infinite square tube of cross section A, filled with 
classical Yang-Mills Higgs plasma at the equilibrium temperature. Two semi-infinite 
regions of definite phase are separated by a domain wall. Because the system has 
translational invariance, for times much greater than any thermalization time in the 
plasma the position x of this wall will diffuse, 

,37) 

with D;j. the diffusion constant for x. Hence, the probability that, starting from x(0) 
at time 0, one will arrive at x{t) at time t, which we call V{x{t),x{0),t), satisfies 

dx{t)Vix{t), x{0),t){x{t) - x{0)Y = \t\D^ . (38) 

Now suppose that we exert a very small force f{x) = —V'{x) on the bubble 
wall. We choose V to be constant in a large neighborhood of the origin, although V 
eventually turns up so that it is bounded from below and goes to infinity at x — ±oo. 
At leading order the hopping probability will be modified by an offset, 

/ dx{t)Vf{x{t),x{0),t){x{t) - x{0)) = {x{t)-x{0)) = Cft = vM)t (39) 
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with C a constant to be determined. Also, the equihbrium probabihty distribution of 
X will be multiplied by a Boltzmann factor exp(— V^/T) = exp{fx/T) ~ 1 + fx/T + 
{fxY/2T^. Starting with a thermal distribution of wall positions at time t = 0, the 
probability that the wall should be at position x{t) = at time t, relative to its 
probability for starting there, is 

1 = j dx{0)Vf(0, x(0),t) i^l + ^ + j = 1 - ^ + , (40) 

and hence 

C = — . (41) 
The pressure on the wall is P — f / A, so the friction coefficient is 

_ P _ P _ 2T 

"^'^'Jc'dJ.- ^^^^ 



In physical units. 



i^.phys = ^^ and Ap,,, = -^, (43) 

so ?7 oc g^T^, as stated above. To compute the bosonic contribution to the friction 
coefficient at leading parametric order it is then only necessary to determine D^A for 
a bubble wall in a classical simulation. 

It is necessary to make a few changes to the above ideas to implement them in 
a realistic simulation. The simulation occurs in a box with symmetric boundary 
conditions; there are two bubble surfaces, and there is a slight mutual interaction 
between them, depending on their separation; for instance, the motion of one wall 
releases latent heat, which slightly changes the temperature and thereby induces a 
response from the other wall; so even if each wall's position can be described by a 
single coordinate, the motion of the two coordinates will be interdependent. The 
interdependences appear in the motion of the difference of the wall positions, but 
translational invariance ensures diffusive motion of the average of the wall coordi- 
nates. The above arguments apply, except that the force / must be the total force 
on both walls, and D should be computed for the average coordinate. The same final 
expression holds, except that A now represents the sum of the areas of the two walls, 
twice the cross section of the box. In terms of the cross-section of the box and the 
diffusion rate of the average coordinate Dg^y, 

5.3 Numerical results 

We determine the wall position as follows. Every few lattice lengths in time, we 
average over the short directions of a long rectangular lattice, producing a c- 
number function of the long direction, this is smoothed with a Gaussian broad enough 



20 



lattirp sizp 


dr 


t (]attirp iinitsi 


V)„ (lattirp iinitsi 


71 


16^ X 192 


8 


1100 


1.21 ±0.24 


.026 ± .005^'^r^ 


322 X 


8 


640 


.25 ±0.05 


.031 ± .006^6r^ 


362 X 144 


6 


674 


.095 ±0.015 


.027 ± .004^6^-4 



Table 1: Data from computing the friction on the bubble wall. Three runs on three 
lattice sizes and at two lattice coarsenesses are within error. 

to make the phases very clearly distinct but well narrower than the wall. For /?l = 8 
and \l = 0.2 we found a = 3 lattice units sufficient, but for larger where the phase 
transition is weaker the smoothing would need to be stronger (and the cross-section 
of the box would need to be wider). We set a threshold at some fraction of the way 
between the average broken phase value and the average symmetric phase value of 
(typically 30%) and, starting at the point of minimum move out in either 
direction and identify the wall as the first point where the threshold is exceeded. This 
works well because, as noted earlier, the fiuctuations in are much smaller in the 
symmetric phase than in the broken phase. 

In Figure |^ we plot the average of the positions of the two walls in a 16 x 16 x 192 
box. The data are taken every 5 lattice lengths for a period of 1100 lattice lengths. 
We extract a diffusion constant from this data by the method discussed in Appendix 
A, which uses the fact that the coefficients of a sine transform of a Brownian process 
are independent and Gaussian, with variance oc l/fc^. The power spectrum is plotted 
in Figure ^ which clearly shows this powerlaw behavior. Results for this run, for 
a 640 lattice length run in a 32 x 32 x 192 box, and for a 674 lattice length run 
in a 36 X 36 X 144 box at = 6, are presented in Table |l|. In each case several 
(order 10) of the lowest frequency Fourier modes were removed, as discussed below, 
and around 100 of the next lowest modes were used; all higher frequency modes, 
which are more polluted by small errors in the wall finding algorithm and which may 
contain physics on time scales too short for the diffusion description to be valid, are 
also excluded. We checked that 50% changes to these cuts did not significantly affect 
the results. The error bars reflect the quality of the flt, see Appendix A, and are 
typically twice the change resulting from a 50% change in the cuts. The runs are 
mutually consistent, which is a good test since they were performed on lattices of 
different sizes and coarsenesses. For comparison, the friction expected from the one 
loop (free scattering) estimate, from transverse W bosons alone, and using the lattice 
value of = 1.54(7, is .069g^T^, larger by a factor of about 2.5. One expected the 
free particle estimate to exaggerate the friction, because scatterings will thermalize 
the particles as they reflect from the wall, reducing the friction; but the amplitude of 
this reduction could not previously be calculated rehably. 

We should mention a serious note of caution in establishing these results and in 
measuring the friction on a bubble wall in a classical simulation. This has to do with 
hydrodynamics. As noted earlier, the liberation of latent heat by a moving bubble 
wall raises the temperature of the plasma locally, and the heat is redistributed by the 
plasma. The temperature at the bubble wall is generically elevated with respect to 
the temperature at the time the bubble nucleated, slowing the wall in a way which 
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cannot strictly be said to be frictional. In the quantum system heat is transported 
by hydrodynamic waves which travel at the speed of sound Vs = l/\/3, and for an 
isolated, spherical bubble, the dissipation of heat into the surrounding medium is 
efficient if the wall velocity is small. 

Latent heat is also liberated in the classical lattice system, but the specifics of 
the hydrodynamics are cutoff dependent and qualitatively different. For instance, in 
physical units the heat capacity of the classical plasma increases with inverse lattice 
spacing as a~^, while the latent heat depends on a~^; so the relation between the 
latent heat and the energy drop between equilibrium and nucleation depends on a. 
Worse, heat is not conveyed on the lattice by bulk hydrodynamic flow; it travels 
diffusively. The reason is that the lattice does not respect momentum conservation; 
while momentum is approximately conserved for interactions between infrared modes, 
many ultraviolet modes obey dispersion relations in which adding momentum acceler- 
ates the mode in the opposite direction of the applied momentum. Hence, scatterings 
cause momentum to be lost into the lattice, and bulk flows tend to come to a stop; 
on large length scales heat dissipates rather than flowing hydrodynamically. We have 
verifled this numerically by heating the lattice system inhomogeneously so as to ex- 
cite the fundamental standing wave; the amplitude of the standing wave decays rather 
than oscillates. 

The problem lies in preventing the hydrodynamic response of the classical fluid 
from contributing to the observed friction. In the case of the diffusion rate mentioned 
above, this becomes a problem for the very low frequency Fourier modes. The latent 
heat liberated by the motion of the wall tends to sit on top of the wall, diffusing 
away only slowly, so there is a very long term "memory effect" inducing negative 
correlations in the wall motion. In Figure ^ we see that the most infrared sine 
transform coefficients do indeed fall off the power law curve; this is why it is necessary 
to drop the most infrared coefficients from the analysis of the diffusion constant. 
Between the low frequency coefficients which must be removed and the high frequency 
coefficients which do not show diffusive behavior there is a sizeable fiducial range, as 
we established by moving the cuts in and out on either side by 50% and finding little 
difference to the determined diffusion rate, as discussed above. The cut for (3^ = 6 
had to be sharper because the heat capacity is smaller compared to the latent heat 
in that case. 

Another way of measuring the infrared contribution to the bubble wall friction 
is to supercool the system in the symmetric phase and produce a small region of 
broken phase by some means, for instance by briefly making mj^Q more negative in 
a narrow region, or by starting with a mixed phase conflguration and applying the 
cooling very abruptly. One then directly measures the velocity at which the bubble 
walls move and sweep up the symmetric phase, averaging over several repitions with 
different initial conditions to determine and improve statistics. This technique has 
the advantage that it does not rely on a linear approximation, which may not be 
very good if the departure from equilibrium is large; but the latent heat poses a 
problem for this technique. For the lattice coarseness and value of Al used in this 
paper, the temperature rise due to the liberation of latent heat is comparable to the 
supercooling, as can be seen from how much the temperature jumps on changing 



22 



phase in the hysteresis plots, see Figure |^. Another problem is that, if the departure 
from linear response really is important, then the frictive pressure may not be linear 
in Vui, and it becomes more difficult to add the effects of the top quark and the W 
boson. 

We end this section by commenting on possible ways around the hydrodynamics 
problem. One way is to work on a much finer lattice, so that the heat capacity is large 
enough to soak up the latent heat without much effect. The obvious disadvantage 
is that this technique is numerically expensive. Alternately, one could add extra 
fundamental representation, massive scalars to the theory to serve as an additional 
heat reservoir. If their masses are chosen large enough then they will have little effect 
on the infrared physics. They would also tend to raise the Debye screening mass, 
which might be desirable. As long as they had no coupling to the Higgs boson, one 
would not expect them to contribute to the friction, except by their interactions with 
the other particles. A final possible solution is to turn on a very weak Langevin noise 
and damping term. The technology for doing this in a gauge invariant way was worked 
out in As a thermalization algorithm the technique is non-optimal because it 
thermalizes the infrared modes much more slowly than the ultraviolet ones; but for 
this application that is a blessing, since one wants to absorb the heat going into the 
ultraviolet excitations without interfering with the correct infrared dynamics. The 
main disadvantage is that the Langevin term makes the evolution canonical, rather 
than microcanonical, so the two phase configuration would be unstable; but in a large 
enough box it would take a long time for the walls to diffuse into each other. We will 
not pursue any of these techniques here. 



6 Chern-Simons number motion on the bubble wall 

Another question which the real time technique can answer is how the system responds 
to a chemical potential or other driving force for the motion of Chern-Simons number 
Ncs- This is directly relevant to the study of baryon number violation because the 
axial anomaly relates the motion of baryon number to the motion of Nqs- In fact the 
classical technique was originally proposed to investigate this question , and since 



then it has been used to investigate Ncs violation in Yang-Mills theory ||2^, in 



the symmetric phase of Yang-Mills Higgs theory |18, M, and in the broken phase of 



Yang-Mills Higgs theory |2^, . It has yet to be applied to perhaps the most relevant 



problem, which is the rate of Ncs motion in the out of equilibrium environment of a 
moving bubble wall during the phase transition. The study of mechanisms for baryon 
number violation at the electroweak phase transition is the study of how the out 
of equilibrium phase transition physics can induce in the infrared bosonic effective 
theory CP violating operators. To convert this information into a baryon number 
abundance we must understand how Ncs responds to these operators. 

Turok and Zadrozny ||^ and McLerran et. al. take the fermions to be in 
equilibrium and integrate them out, and then try to investigate how the resulting 
CP violating operators will influence the out of equilibrium decay of gauge fleld 
configurations as they are swept onto the wall; but for lack of quantitative tools 



23 



they were forced to make only qualitative or parametric estimates (see also 
More recently work has focused on how the transport of fermions can carry CP 



violation away from the bubble wall [Q, (^Of, into the symmetric phase. One of the 



main motivations for examining the transport mechanism was the worry that baryon 
number violation on the wall would be strongly suppressed except for at the very 
leading edge 0, ||4^, by the turning on of the sphaleron mass. Thus a chemical 
potential inside the wall would be ineffective. But it remains to be shown whether 
this picture, based on viewing the baryon number violation as a completely local, 
quasiequilibrium process, is correct, or whether decaying gauge field configurations 
swept up onto a bubble wall can still respond to CP violating effects on the wall to 
produce Ncs change. In this chapter we will focus on this question. 

To answer it we need to apply a chemical potential for Nqs on the lattice, as 
otherwise we cannot create a CP violating biasj]. The technology for doing this has 
recently been worked out for Yang-Mills theory [^, and it is straightforward to 
extend it to Yang-Mills Higgs theory. We briefly review the method in Appendix A, 
and extend it to Yang-Mills Higgs theory. For a complete exposition see 

The only problem in applying a chemical potential for Ncs in Yang-Mills Higgs 
theory is that the rate of Nqs niotion (or the rate of Ncs diffusion when no chemical 
potential is applied) is slightly contaminated by ultraviolet lattice artefacts; in the 
broken phase these effects contribute about 1/10 of the symmetric phase rate to the 
response. We give evidence that this rate is an ultraviolet lattice artefact in Appendix 

A. The consequence of this problem is that we must observe a rate of Ncs motion well 
above 1/10 the symmetric phase response rate before we know that it corresponds to 
genuine infrared processes and not to spurious ultraviolet physics. 

With this in mind, we set out to investigate how Ncs interpolates between the 
symmetric phase value and the broken phase value across a bubble wall, when there 
is a spatially uniform chemical potential for Ncs in both phases. We do this first in 
equilibrium. We take a mixed phase configuration with about equal volumes of broken 
and symmetric phases, in a 16 x 16 x 192 box at the equilibrium temperature and 
the values of = 0.2, = —.3223 generally used in this paper. We compute Ncs 
at each lattice point and bin it according to its distance from the bubble wall. This 
means that we find the bubble wall surfaces by the algorithm described in Appendix 

B, and for each point in the plasma we find the minimum vertical distance to a bubble 
wall, and add Ncs to a bin corresponding to that distance, considered positive if the 
point is on the broken phase side of the wall and negative if it is on the symmetric 
phase side. Binning this way accounts for the fact that the wall is an uneven surface. 
We also bin <l>^$, unsmoothed, using the same rule; this gives a bubble wall profile. 

The results for the equilibrium case with a space independent chemical potential 
are presented in Figure The Ncs results have been smoothed with a Gaussian 

^Another possibility is to put CP violating high dimension operators, such as ^^^E ■ B, in the 
action; although such terms cannot bias Ncs in equilibrium, the shift in <i>^<l> during the phase 
transition will bias topology change. However, implementing such operators makes evolving the 
system much more complicated, because the update rule becomes nondiagonal in the natural basis of 
degrees of freedom; also any nonrenormalizeable operator has much more influence on the ultraviolet 
behavior of the lattice system than on the infrared, which is potentially dangerous. 
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envelope of width a = 3 lattice units to eliminate white noise fluctuations, while the 
plot is completely unsmoothed; its smoothness arises from averaging over a long 
time (approximately 4000 lattice lengths) and over the area of the bubble walls. What 
we see in these plots is that, in equilibrium, the response to a chemical potential falls 
off abruptly just inside the bubble wall, and sphaleron events are suppressed in most 
of the interior of the wall, as well as behind it. 

Next we consider the out of equilibrium case. We produce a series of initial condi- 
tions by evolving a mixed phase, equilibrium conflguration through a series of short 
{t = 10 lattice unit) Hamiltonian trajectories. The starting conflguration contains 
mainly symmetric phase, so the broken phase can expand for some distance before 
the walls meet. For each initial condition, we evolve it using the "local thermostat" 
discussed in Section |] to drive the temperature to a point midway between equilib- 
rium and the spinodal point; we also turn on a global chemical potential for Nqs- 
After waiting about t = 20 lattice units for the system to achieve a steady state, we 
begin to record Ncs and binning according to proximity to the bubble wall. We 
discontinue the evolution well before the bubble walls meet. 

The wall shape for the out of equilibrium case, and the motion of Ncs in this 
case, with a constant chemical potential, are presented in Figure ^. The conclusion 
is similar to the case of the wall at rest; baryon number shuts off some short distance 
into the bubble wall. 

We also investigate the case where the chemical potential exists only on the bub- 
ble wall. We could do this by adding a nonrenormalizeable operator ■ B to the 
action, which closely resembles an operator which would arise in the two Higgs dou- 
blet model. However there are serious technical difficulties with adding such a term; 
it significantly modifies the ultraviolet dynamics, and it renders the update rule non- 
diagonal. Instead we will "mock up" the effect of this term. First, we measure the 
averaged dependence of <I>^$ on distance from the bubble wall, which is shown in Fig- 
ure p. Now for each point in the plasma, we find the vertical distance to the nearest 
bubble wall, and apply a chemical potential at that point which is proportional to 
the derivative of the wall profile in Figure |^ at that value. The chemical potential is 
then only nonzero on the wall, and in a way which mimics d{^^^)/dz, which would 
be oc d{^^^)/dt for steady motion of the bubble wall. 

We have performed a series of runs with a chemical potential which is only nonzero 
"on the wall" . The bubble wall shape and the average of Nqs for each run are plotted 



in Figure |T0[ The figure shows that Ncs is generated in the interior of the bubble 
wall, where the chemical potential is being applied, but it is destroyed in a region 
immediately to either side. This is not a fluctuation due to insufficient statistics; 
if the data set is split in four, each quarter shows the same morphology. The rate 
of baryon number generation within the hump itself is, in the dimensionless units 
discussed in Appendix A, 

iPL7ryk.r.,Ncsd'^dt (45) 
/ ^LfJ'ix, t)d^xdt 

However, when one integrates over the full simulation volume, including the regions 
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to either side of the bubble wall where Ncs is negative, one finds 



((3lt^)U Ncsd^xdt , , 

= 0.091 ±0.025. (46) 

The error bar on the latter number is estimated assuming the error is from diffusion 
of Ncs in the symmetric phase volume in the simulation, and is consistent with the 
error bar from statistics between runs, the final result is larger than the expected 
"UV artefact rate" of kuv artefact — 0.055 (see Appendix A), so one might conclude 
that there is a net infrared generation of Nqs from the chemical potential; however 
the statistical significance is not very good, and the generation is about an order of 
magnitude less efficient than in the symmetric phase. 

Why is there a bump with a dip on either side, and how can the system respond to 
a chemical potential on the wall? Recall why a chemical potential should not induce 
baryon number violation in the broken phase. If a chemical potential is turned on 
briefly, it will push the gauge flelds in the direction which generates baryon number, 
and, briefly, Nqs will be generated. However, the flelds will respond elastically to 
this impulse-after oscillating in the direction of greater Ncs they will "bounce back" 
to the previous value of Ncs- The temporal response to an impulse of chemical 
potential will then be a forward surge in Ncs followed by a resotration to the original 
value. (We have tested this statement by tracking Ncs after such an impulse.) Now 
the gauge flelds also propagate, so if a chemical potential is applied in a spatially 
nonuniform way, the relaxation to the former value of Ncs may take place not where 
the chemical potential was applied, but nearby. This would explain the generation of 
Ncs on the wall and the destruction of Ncs immediately to each side. In the case of 
the moving wall, some of the "excited" gauge fleld modes relax inside the symmetric 
phase, because they propagate there. Once in the symmetric phase, their behavior 
can be different-without a Higgs condensate present, the most infrared modes may no 
longer be oscillatory and will not relax as much as they would have, so the deficit in 
Ncs on the sides of the wall need not be as large as the production on the wall. The 
infrared gauge fields themselves can serve to transport the CP violation on the bubble 
wall to the symmetric phase, where it can lead to baryon number nonconservation. 
However, our data show that this is quite an inefficient mechanism. 

In a previous draft of this paper we found a considerably larger effect from chemical 
potential on the bubble wall. However, in the evolutions used for those results we 
were insufficiently careful to stop runs well before the bubble walls met; sometimes 
they would begin to collide, causing the bubble wall finding algorithm to go awry, 
and possibly leading to the application of chemical potential inside the symmetric 
phase rather than on the bubble wall surface. For the runs presented here we were 
more careful; the new runs also represent about 3 times the cumulative evolution time 
(totahng t — 14000 lattice units). 
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7 Conclusion 



We have reviewed the fact that the thermodynamics of quantum Yang-Mills Higgs 
theory coincides, in the approximation of dimensional reduction, with the thermody- 
namics of the classical theory, and shown how this can turn the real time evolution of 
the classical theory into a powerful microcanonical Monte- Carlo algorithm. We used 
the algorithm, for Higgs potential parameters corresponding to a tree level vacuum 
Higgs mass of ~ 50GeV, to investigate the phase diagram, including both metastable 
branches, and to find the equilibrium temperature. We also developed and applied 
technology for finding the equilibrium surface tension by microcanonical techniques. 

Further, we have strengthened the argument that the dynamics, as well as the 
thermodynamics, of the infrared can be approximated with classical field theory; the 
errors should be parametrically suppressed (order aw) for phenomena which are truly 
infrared dominated. This is a consequence of the weakly coupled nature of the theory; 
the physics can only become nonperturbative by becoming classical. We were thus 
able to compute to leading parametric order the friction from bosons on a moving 
bubble wall in the near equilibrium limit, and to investigate the physics of baryon 
number violation out of equilibrium in the presence of a moving bubble wall. We have 
found that the friction from infrared bosons is smaller than that from transverse W 
bosons in the free scattering limit by a factor of 2 - 3, and baryon number violating 
processes proceed some distance into the bubble wall but stop well short of the broken 
phase. 

The classical method opens all thermodynamic and dynamical properties which 
are dominated by the infrared bosonic sector to reasonably accurate calculation. To 
pursue the thermodynamics with good accuracy it will be necessary either to apply 
a great deal of computer time, or to make a thorough investigation of linear in a 
corrections. To improve the accuracy and reliability of the dynamical calculations 
it will be useful in addition to find some way to properly represent hydrodynamic 
effects on the lattice. It may also be possible to integrate out the field, so the 
large Dcbyc mass limit is obtained without necessitating very fine lattices. 

We are hopeful that such improvements in accuracy can be achieved, making it 
practical to use the techniques developed here at a wide range of Higgs masses and 
for interesting extensions of the Minimal Standard Model, including a second light 
Higgs doublet and a fight stop squark. 
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A Ncs in Yang- Mills Higgs theory 



Two quantities can be used to characterize baryon number violation in classical Yang- 
Mills Higgs theory; the diffusion rate of Ncs, and Ncs when there is a chemical 
potential for A^^^ added to the Hamiltonian. In the large volume, long time limit 
these quantities can be characterized by two rates. 



lim 

t^oo 



{{Ncsjt) - NcsjO)?) 
Vt 



and 



r„ = lim 



NcsT 



(47) 



2Tfj_, derived in 



|2 



and further 
if Ncs rises 



which satisfy a fluctuation dissipation relation, T 

discussed in |]S3|, It is also interesting to know Ncs ior fi ^ T 
much faster than linearly in fi it indicates that the motion of Ncs is obstructed by a 
substantial free energy barrier; for Yang-Mills theory Ncs turns out to be very linear 



so there is no free energy barrier [21 



The technology for computing N, 



cs as a function of time in Yang-Mills Higgs 

1]; we will give 



theory, in the absence of a chemical potential, was developed in |]T9 
only the most cursory review here. Instead this appendix will concentrate on how 
best to extract a diffusion constant from the diffusive trajectory of Ncs, and then 
on how to extend the chemical potential method developed in [^] to the case of 
Yang- Mills Higgs theory. We will then investigate Ncs motion in the symmetric and 
broken phases, using both techniques, and present three pieces of evidence that the 
rate we establish in the broken phase arises from spurious ultraviolet effects. 



A.l Finding a diffusion rate 

Suppose we know some quantity 2; as a function of time t between an initial time 
and a final time tf. We believe that it should exhibit Brownian motion, and we want 
to determine the diffusion constant F, defined as {{z(ti) — z(t2)Y) = T\ti — 12\- One 
method, used in [^, is to compute 

lag(ti) = r '\z{t + ti)- z{t)fdt (48) 



for many values of ti and to fit the result to a straight line. While this method works, 
the fitting is complicated by the large correlations between lag(t) for different values 
oft. 

We advocate an alternate method in which one sine transforms z and makes a 
likelihood analysis of the transform coefficients. The analysis is simpler because the 
sine transform coefficients are independent; and the likelihood analysis is probably 
optimal in the sense of statistical power. 

We begin by redefining z{t) as z{t) — -2(0); then 

{z{ti)z{t2)) = Tm.m{t,,t2). (49) 

Defining the sine transform coefficient 

*f dt , ^ ( mxt\ , , 

-z(t)sinf— j n = 1,3,5,7... (50) 
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we quickly find 




The sine transform coefficients are independent, with variance ^Ttf/ (nvr)^. They will 
also be Gaussian distributed. We can see that the sines used form a complete set by 
extending z as an odd function about the origin and an even function about t/. 

In the case of interest, z will only be known at discrete values of t, and the in- 
tegrals above represent the corresponding sums, z will also often be contaminated 
by white noise-this was proven in the case of Yang-Mills theory by Ambjorn and 
Krasnitz, by considering the Abelian approximation to the ultraviolet modes of the 
theory [pO| -which will introduce a constant times 6mn into Eq. (|5^) . The value of 
this constant is uninteresting, but its presence means that the ultraviolet coefficients 
do not carry information about F, which can only be determined with finite preci- 
sion. Any other non-Brownian correction to the ultraviolet physics, as we expect for 
instance in extracting the diffusion constant for the bubble wall motion considered in 
Section ^, also makes the ultraviolet data useless. 

The remaining question is how, given a set of Gaussian distributed coefficients 5„ 
, n = 1, 3, 5 . . ., satisfying 

C^'n) = -,+B, (53) 

to extract A and its errorbars. If we assume a uniform prior for the values of A and 
B, then Bayes' theorem gives the likelihood of values A and B as 

ViA,B)oc n ^2iP<l2l/^. .?. = 4 + B. (54) 

n=i,3,... V27ra„ 

The best value of {A, B) is that which maximizes this likelihood function. If the 
likelihood function is sharply peaked and there is a neighborhood of the maximum 
which contains almost all of the probability, and in which InP is well approximated 
by a quadratic form, then the likelihood will be approximately Gaussian distributed 
with an error matrix which is the inverse of the quadratic form. If there is enough 
data to make a strong determination of A and B then this is generally the case. This 
is how we determine best fits and errors for Brownian processes in this paper, and it 
will allow us to find the Ncs diffusion rate in either electroweak phase. 



A. 2 Chemical potential method in Yang-Mills Higgs theory 

It is useful to be able to measure the response of A'",^^ to a chemical potential. This 
idea was first investigated for Yang-Mills Higgs theory by Ambjorn et. al. [|19|, who 



encountered technical difficulties; these were overcome for Yang-Mills theory in 
Here we briefly extend the method developed there to Yang-Mills Higgs theory. 
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The first step is to define an (adjoint valued vector) magnetic field B; the appro- 
priate choice is 



BK^) = ^ E ^Tr - zt'^Uu + 1 E ^Tr - zU,{x)T'^Uj{x)Un (55) 

where each sum is over 4 plaquettes orthogonal to the x, i link, [/□ is the product of 
links around that plaquette beginning and ending on the link, and the two sums are 
for the basepoint and endpoint of the link. The Ui{x) and Uj{x) in the second term 
parallel transport [/□ to the basepoint of the link. Ncs is then defined as 

27T^Ncs = Y.E-Bf = E-B (56) 

where we have displayed the definition of inner product for adjoint vector fields. These 
two definitions are already sufficient for tracking the (diffusive) motion of Ncs in the 
absence of a chemical potential. 

Next we see how to apply a chemical potential for Ncs- First consider Yang- Mills 
theory. On the lattice it is not necessarily true that 

= V ■ B{x) = Y.iBi{x) - U}{x - i)Bi{x - i)Ui{x - i)) (57) 

i 

which means that if, for some reason, the Gauss constraint were not actually obeyed, 
Ncs would depend on the size of the violation. To understand this point, suppose 
that we knew an orthonormal basis for adjoint vector fields which partitions into basis 
elements which are linear combinations of the constraints (which can be constructed 
by taking the gradients of adjoint valued scalar fields) and basis elements with zero 
divergence. Call the projection of the electric field into the former subset E"^ (c for 
constraint) and the projection into the latter subset E*. The Gauss constraint is the 
condition E'^ = 0, and the E* are the dynamical degrees of freedom. Similarly, B 
partitions into B^ and B*, and the above point is that B^ ^ 0, so 

2Tr^Ncs = E^ ■ B" + E* ■ B* (58) 

would change value if for some reason one orthogonally departed from the constraint 
manifold. 

This is the source of problems for the implementation of a chemical potential for 
Ncs, which modifies the E equation of motion by 

Since i?^ 7^ this moves E'^ away from zero, which in turn changes the value of Ncs, 
producing wrong answers. What has happened is that we have added a new term to 
the Hamiltonian, with a part linear in E'^, so the constraint is no longer first class. 
The evolution then departs from the constraint manifold. But if we redefine Ncs as 

27r^Ncs = E*-B\ (60) 



30 



which coincides with the previous definition on the constraint manifold, then the 
constraint will again be first class and the evolution will not violate Gauss' law. The 
modification to the E field equation of motion becomes 

SE* = -J^B*, 6E^ = 0, (61) 

which preserves the Gauss constraints. And since we never leave the constraint sur- 
face, we can still compute Ncs using E ■ B. 

The easiest way to implement this addition to the equation of motion is to apply 
Eq. ( ^9] ) but to then remove the contribution to E^ by orthogonally projecting to 
the constraint surface. An exact orthogonal projection is numerically expensive; a 
very accurate approximate algorithm, which is exactly orthogonal but does not quite 
complete the projection, was developed in (which is where the interested reader 
should go for more thorough details). It is also shown there that this technique gives 
a value for equal to rd/2, in accord with the fiuctuation dissipation relation, a 
good check. 

In the case of Yang-Mills Higgs theory, it is no longer true that E^ = 0; instead 
the Gauss constraints stipulate that 



0. 



Wx, a - Re (u\x)iT''^x)) + ^ [e^{x) - {Uj{x - i)Ei{x - i)Ui{x - i)Y 

(62) 

To analyze the application of a chemical potential it is convenient to consider an 
orthonormal basis for the momenta E, 11 which partitions into 4 subcatagories: the 
E*] the radial components of 11; an orthonormal basis of the linear combinations 
of degrees of freedom which are forced zero by the Gauss constraints, P^; and an 
orthonormal basis of the remaining degrees of freedom, P*. and P* are mixtures 
of the E'^ and Higgs momentum degrees of freedom. The definition of Ncs can be 
written as 

2ti'^Ncs = B ■ E = B* ■ E* + ■ P* + B" ■ P" , (63) 

where the dot products involving the P are between B and the E field components of 
the P. As above, applying a chemical potential using this definition of Ncs will excite 
violations of the Gauss constraints; again the solution is to say that the definition of 
Ncs is 

27r^ Ncs = B* ■ E* + B"" ■ P* , (64) 

which is equivalent on the constraint manifold and therefore does not require a change 
in the program code which computes it. All that changes is that the momenta should 
be orthogonally projected to the constraint manifold each update. The algorithm for 
this orthogonal projection is the same as the algorithm proposed for use in thermaliz- 



ing Yang-Mills Higgs theory in ||21| and used in this paper for the canonical ensemble 
algorithm. 

A. 3 Results for Ncs motion in Yang-Mills Higgs theory 

Now we will apply these tools to investigate the motion of Ncs in Yang-Mills Higgs 
theory in each phase. Data for the diffusion rate in Yang-Mills and Yang-Mills Higgs 
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Theory 






Phase 


t 




"expected errorbar" 


YM 


6 


N.A. 


N.A. 


1000 


1.05 ± .11 


±.18 


YMH 


~ 8 


-.3223 


symmetric 


4000 


0.85 ± .07 


±.14 


YMH 


~ 8 


-.3223 


broken 


8000 


0.111 ± .009 


±.037 


YMH 


8 


-.50 


deep broken 


3000 


0.059 ± .008 


±.04 



Table 2: Diffusion rates of Nqs in Yang-Mills and Yang-Mills Higgs theory. In the 
last column the bare Higgs mass squared was rriHo = —0.5 in lattice units, placing 
the system deep in the broken phase. The "expected errorbar" is the errorbar which 
would occur if the diffusion were made out of Poisson distributed integer steps. 



theory are given in Table 0. As has become standard in the literature, the rate is 
expressed in terms of 

= r£;(7r/3i)^ lattice units or rrf(atf,T)~'' physical units . (65) 

All data are for 16^ lattices, which should be large enough to achieve the infinite 
volume limit for the values of (3^ used The Yang-Mills Higgs data are all for 

(bare) = 0.20. To illustrate the data analysis technique. Figure ^ shows a 6000 
lattice unit section of the evolution of the Yang-Mills Higgs system in the broken 



phase, and Figure ^ gives the coefficients of the sine transform, clearly showing the 
A/n^ + B behavior as well as the wide scatter associated with the log of the square 
of a Gaussian quantity. 

There are three pieces of evidence which make it difficult to believe that the ob- 
served diffusion constant in the broken phase represents genuine infrared, topology 
changing physics. The first is rather subtle; it involves error bars. A series of integer 
steps in Ncs would not appear as perfect diffusive motion, and this would be re- 
flected in the error bars in the extracted diffusion constant; if the steps were Poisson 
distributed, as seems reasonable, then N steps would give error bars in of F^/ a//V. 
The number of steps, based on F^, is TaVt, and the error estimate based on this 
reasoning is shown in the last column of the table. In the broken phase the actual 
error bars are much too small, suggesting that the diffusive process is somehow much 
smoother than expected. This would happen if it was the accumulation of a large set 
of small, ultraviolet shifts in Ncs- 

Ambjorn and Krasnitz have shown that, at leading order, the ultraviolet theory 
behaves like an abelian theory, and shifts in Ncs are elastic, and are restored a 



moment later pO[. But they also show that such shifts occur at a rate which diverges 



as a — > 0, so if at next to leading order such shifts are not perfectly restored, it 
could still contribute an effect which would not vanish as a — > 0. This motivates the 
possibility of a spurious signal which remains constant in the small a limit, which 
would neatly explain the behavior of the diffusion rate observed above. 

For a further piece of evidence, we applied a very large chemical potential, fi = 
6/Pl, to a 16^ lattice of Yang-Mills Higgs plasma in the broken phase just below 
the phase transition temperature. Over a time period of t = 4000, Ncs changed by 
11.4, which would correspond to Kd = .088 ± .026 if Ncs rises linearly with /i out to 
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6//3l. Since this agrees with the measured value of Kd, we conclude that the response 
is indeed linear. However, if there were a substantial free energy barrier we would 
expect the rate to rise with sinh(2/3i:,/i), which it clearly does not. Since symmetry 
is broken strongly enough at the value of we are considering that there should 
be a substantial free energy barrier, this implies that the effect we observe is not 
the system jumping that barrier, but some other process, presumably the ultraviolet 
effect suggested above. This evidence is in the same spirit as the first piece. 

For a final, strong piece of evidence, consider the last line in the table. Here we 
evolved a volume of Yang-Mills Higgs fluid with a large negative Higgs mass squared, 
so it was very deep in the broken phase; in fact we measured = 98, over 5 times 

its value in the broken phase just after the phase transition for the parameter values 
we have concentrated on in most of this paper. This corresponds, in physical units, 
and after subtracting off the ultraviolet 1/a contribution, to = 4.7 gT. The tree 
level Sphaleron energy for this value of is 4.7 x 4:TtB{X/ g'^)T ~ lOOT ||5^. For such 



a large value, the semiclassical Sphaleron approximation should be very reasonable, 
and the diffusion rate should have an immense exponential suppression, by a factor on 
order exp(— 100) [^, yet the measured rate is only mildly lower than that barely 
inside the broken phase. This diffusion of Ncs must be a lattice artefact, unless our 
understanding of baryon number violation deep in the broken phase is completely 
wrong. The slight decline in the rate relative to the rate "just" inside the broken 
phase probably comes about because is larger even than the lattice inverse spacing, 
so even the ultraviolet lattice modes have their populations mildly suppressed. 

We should point out that, while the error bars for diffusion in Yang-Mills theory 
and in the symmetric phase of Yang-Mills Higgs theory are also smaller than expected, 
they are not as much smaller, and the argument that Ncs should diffuse in integer 
steps is invalid in the unbroken phase. Further, it has been demonstrated that the 
diffusion rate falls off for small volumes, which shows that the dominant contribution 
is indeed infrared and should reflect real physics pOf. It is likely that a small part 



~ 10% of the measured rate in the symmetric phase and in Yang-Mills theory arises 
from ultraviolet artefacts; in fact the results here could be viewed as a calibration 
of that contribution, since the ultraviolet behavior of lattice Yang-Mills Higgs theory 
should be about the same in each phase. 



B Bubble wall surface tension 

In this appendix we present the details of the calculation of the bubble wall surface 
tension. The determination has two parts: first we identify the bubble wall surfaces, 
then we show how to extract the surface tension from the shape of the surface (aver- 
aging over a large sample of surfaces). The reader can skip the first subsection and 
just take the wall surface position to be known if they are uninterested in numerical 
details. 
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B.l Finding the surface 

Suppose we have the values of $ and the connections U at some point in time and we 
want to determine where the phase boundary is. The problem is that $ contains a lot 
of ultraviolet fluctuations, so a simple threshold definition of the two phases does not 
work. We will need to coarse grain or smooth the fields. Define the once smoothed 
Higgs field $i as 

$i(x) = ]^x) + ^ V \Ui{x)^x + + U}{x - i)^x - i)] . (66) 

The sum just averages over nearest neighbors, parallel transporting along the shortest 
path. $n+i is defined by applying the same averaging process to $„. For large n, 

is approximately $ averaged over a Gaussian envelope of variance 3n/4. The 
averaging very quickly removes ultraviolet fluctuations. It also very slowly degrades 
the condensate, because the averaging includes an averaging over different parallel 
transport paths, and the trace of any Wilson loop is less than Tr/. After several 
smoothings (we use $§, but for a finer lattice or a larger value of we would need to 
use more) the phases become much more distinct, although one still cannot distinguish 
them by setting a threshold. 

To work with a real c number quantity we take $]j$n- Label the three directions 
Xi,X2i and X3, with 2:3 the long direction in the lattice. For each column (fixed xi 
and X2) we smooth along the column with a Gaussian envelope of width large enough 
to make the two phases distinguishable by a threshold, but smaller than the typical 
wall thickness (which for our value of and /?l will turn out to be on order 15). In 
practice we find a width ~ 7 is sufficient. (If it is impossible to smooth enough that the 
two phases can always be distinguished by a threshold without making the smoothing 
length larger than the wall thickness, then one should have used more gauge invariant 
smoothings earlier.) Call the function we get by doing this smoothing /(a;i, 0:2, 2:3). 
The only problem left is to choose a threshold by which to define the phases. We do 
this as follows. First, for each column we find the minimum and maximum values 
of /; average these over the columns; call them fmin and fmax- Wc assume that any 
point with / more than halfway from /^m fo fmax is in the broken phase and any 
point less than 0.2 of the way is in the symmetric phase, and from this definition we 
find the average value of / in each phase. Then we set a threshold some percentage 
of the way between the average symmetric phase value of / and the average broken 
phase value of /; at equilibrium the choice is 30%, reflecting that fluctuations in / 
are much larger in the broken phase; but for a heavily supercooled system we set the 
threshold higher, 40%, because the symmetric phase fluctuations become larger. 

With the threshold set, we deflne the wall height function z{xi,X2) by starting at 
the value of xz in the xi,X2 column where / is minimum and incrementing X3 until 
the flrst point is reached where the threshold is exceeded; z{xi,X2) is chosen as the 
value of ^3 where this happens. The other surface is found by moving backward. 

The definition of the surface described here is imperfect; it assumes that the 
surface is single valued (never bending more than 90° from horizontal) and it smooths 
the surface at a length scale on order the Gaussian envelope radius used to define the 
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smoothed Higgs field $„. However, the surface is not really well defined on length 
scales much shorter than the wall thickness, and we will be most interested in infrared 
scales on which the surface should look relatively fiat, so neither problem will be of 
consequence. 



B.2 Determining the surface tension 

Take the height function z{xi,X2) to be known; how can one determine the surface 
tension from it? 

The free energy of a surface with surface tension a and height function z is the 
area of the surface times the surface tension, 

F = aJ (i^x^l + {Vzf . (67) 

At finite temperature, surface waves will be excited and the partition function for the 
appearance of the surface will be 

Z = y r>^exp(-/3F), (68) 

which is the partition function of a simple two dimensional Euclidean field theory. In 
the infrared we can expand F as 

F^ajd'x{l + liVzf - liiVzff + ...). (69) 

The first term contributes a z independent constant and the remaining terms describe 
a free massless field theory plus nonrenormalizeable operators. These operators will 
cause some renormalization of a between an ultraviolet scale and the deep infrared, 
but because they are nonrenormalizeable a will have an unambiguous deep infrared 
limit, which is the value of interest in any calculation where the surface tension is 
a useful concept, concentrating on the infrared and remembering that there may 
be corrections for ultraviolet modes, and taking the area over which the surface is 
stretched to be an L x L square, 

/3i^^^////rf'^(l + ^(V^r) • (70) 
This free field theory is trivially solved by Fourier transformation. Defining 

d^x 



we find 



and by equipartition 



!a = /^^We"«"" (71) 



^F=Z'-^ (72) 



<^->-si?- P3) 
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This expression should be obeyed in the hmit of large L/n; for smaller values 
of L/n it will have corrections because of the nonrenormalizeable operators, and in 
the numerical calculation it will have corrections in the ultraviolet because of the 
properties of the numerically determined bubble wall surface mentioned earlier. To 
compute the surface tension it will then be necessary not only to average the z„ 
over a large sample of walls, but to perform an extrapolation to zero n. Because 
the corrections to the free field theory are nonrenormalizeable, the corrections to Eq. 
([73|) will appear as powers of (ra/L)^; there will be no logarithmic corrections. As 
discussed in the text we performed the large L/n extrapolation by fitting 'n?'{z^) to 
ylexp(— Sn^). We do not know of a good physical motivation for this form except 
that it resembles what a Gaussian smoothing of the wall surface would do to the 5„. 
In practice we could equally well fit only the first few points to a straight line; the 
result is about the same. 
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Tc El 

Temperature Energy 

Figure 1: Schematic depiction of the dependence of an order parameter on temper- 
ature and energy, inchiding metastabihty hues and the equihbrium coexistence line. 
Though equihbrium phase coexistence only occurs precisely at T^, it happens for a 
range of energies between Ei and E2. Phase mixture is stable for a fixed energy 
system if the energy is in this range. 
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Figure 2: Lattice ^t^/j-s ^ ^ function of temperature (I/Pl) for a 30^ lattice; at 
left, without smoothing, and at right, with one nearest neighbor smoothing. The 
system begins in the symmetric phase at high temperature and is cooled through 
the phase transition, and is then heated back to the original temperature. The sys- 
tem supercools during the cooling and superheats during the heating, resulting in a 
hysteresis loop. The coohng represents a total elapsed time of 1280 lattice units; the 
heating is twice as long, with longer bins, to suppress fluctuations. The fluctuations in 
$t<|)/2"2 g^j,g much larger in the broken than the symmetric phase, except near where 
the symmetric phase becomes unstable, where the fluctuations become larger. Both 
the jump in the order parameter and the fluctuations are almost identical before and 
after smoothing, so both are infrared phenomena; but the symmetric phase value of 
is primarily ultraviolet, and is suppressed by the smoothing. 
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Figure 3: Average of the trace of some Wilson loops for the same run. The normal- 
ization is such that the trace is 1 in vacuum. Fluctuations in the Wilson loops are 
about the same size in each phase, and they follow quite closely the fluctuations of 
For the 9x9 loop the symmetric phase Wilson loops are almost zero, showing 
no long range order; they are zero within tight statatistical error for 13 x 13 loops, 
and here the broken phase values are nearly zero as well. 
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Figure 4: Six times smoothed Higgs field, averaged over the short directions of a 
16 X 16 X 192 lattice at the equihbrium temperature T ~ 0.1241, plotted as a function 
of position in the box. The two phases and the phase boundaries are clear, although 
the broken phase has quite sizeable fluctuations in the order parameter. 
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Figure 5: Bubble wall surface power spectrum: {z^'n?) plotted against n^. The upper 
data are for a 32 x 32 cross section box at /3l = 8 and the lower data are for a 36 x 36 
cross section box at [3l = 6. The fit functions are exponentials, used to extrapolate 
the data to the infrared {n — 0) limit. 
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Figure 7: Power spectrum for the previous plot, obtained by sine transformation. 
The wide scatter is expected for the log of the square of a Gaussian distributed 
quantity. If the motion is Brownian then the points should be uncorrelated and 
should depend on n as 1/n?. 
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Figure 8: Wall shape in equilibrium (solid line). The vertical axis is in lattice 

units, equal to 4(0^/g(^T^) in the continuum. The symmetric phase value has been 
subtracted. The horizontal axis is in lattice units, which equal 4/(/3l5'^T) ~ 1.2/T in 
physical units; the zeropoint of the axis is arbitrary. The dashed line shows dNcs/dt 
in reponse to a constant chemical potential, as a function of position relative to this 
wall. The vertical scale for Ncs is arbitrary. The rate falls off sharply inside the wall. 
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Figure 9: Wall shape out of equilibrium, below the phase transition temperature 
(solid line). The wall is in this case moving with a speed of f ~ 0.3. The dashed line 
shows dNcs/dt in reponse to a constant chemical potential, as a function of position 
relative to this wall. Again the rate falls off sharply inside the wall. Axes as in 
previous Figure. 



47 




60 80 100 120 140 

position(lattice units) 



Figure 10: Bubble wall shape and Nqs when the chemical potential for Ncs was 
proportional to the gradient of the wall. The vertical axis for Ncs is arbitrary. There 
is a spike in Ncs on the wall, where the chemical potential was applied, and a pit on 
either side of the spike. 
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Figure 11: Chern-Simons number diffusing in the broken electroweak phase in a 
16^ box at — 8. The diffusion does not resemble sudden discrete jumps between 
integer spaced plateaus. 
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Figure 12: Sine transform of the previous data. The spectrum is in excellent agree- 
ment with white noise (constant power) plus a Brownian signal (power oc n~^). As 
discussed in the text, the Brownian motion is probably a lattice artefact. 
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